function [diff, eq] = find_stationary_dist(eq, param,glob,options)
    
    Q           = eq.Q;

    % Distribution of incumbent firms
    L               = eq.L;

    % Distribution of entrant firms
    G = eq.G;
    
    
    % which entrant firms choose to enter?
    
    if (sum(G) < 1e-10)
        L = 0*L;
    else
        for itL = (1:options.itermaxL)
            Lnew    = Q'*((1-param.gamma).*(1-eq.p_exit).*L) + param.M*Q'*G;  
            dL      = norm(Lnew-L)/norm(L);
            if (dL<options.tolL),break,end

            L       = Lnew;
        end
    end

    % L is the distribution over all states
    % Lp is the distribution over customer bases

    Lp = nan(size(glob.bgridf));
    for i = 1:length(glob.bgridf)
        Lp(i) = sum(L(find(abs(glob.sf(:,1) - glob.bgridf(i)) < 1e-6)));
    end

    eq.L  = L;
    eq.Lp = Lp;
    
    q = eq.yf/param.C;
    diff = (sum(eq.L.*upsilonp(q,param).*q))^(-1) - param.D; 

    eq.Dout = (sum(eq.L.*upsilonp(q,param).*q))^(-1);
    eq.mass = sum(param.M*G);
    
end